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We compare dynamics and waveforms from binary neutron star coalescence as computed by 
new long-term (~ 10 orbits) numerical relativity simulations and by the tidal effective-one-body 
(EOB) model including analytical tidal corrections up to second post-Newtonian order (2PN). The 
current analytic knowledge encoded in the tidal EOB model is found to be sufficient to reproduce the 
numerical data up to contact and within their uncertainties. Remarkably, no calibration of any tidal 
EOB free parameters is required, beside those already fitted to binary black holes data. The inclusion 
of 2PN tidal corrections minimizes the differences with the numerical data, but it is not possible to 
significantly distinguish them from the leading-order tidal contribution. The presence of a relevant 
amplification of tidal effects is likely to be excluded, although it can appear as a consequence of 
numerical inaccuracies. We conclude that the tidally-completed effective-one-body model provides 
nowadays the most advanced and accurate tool for modelling gravitational waveforms from binary 
neutron star inspiral up to contact. This work also points out the importance of extensive tests to 
assess the uncertainties of the numerical data, and the potential need of new numerical strategies 
to perform accurate simulations. 



I. INTRODUCTION 

Gravitational waves (GWs) emitted by binary neutron 
star (BNS) inspiral and coalescence will be detectable 
by advanced LIGO- VIRGO detectors. The tidal signa- 
ture in such waves is (mainly) proportional to the tidal 
polarizability parameter fj,2 that yields the ratio between 
the tidally induced quadrupole moment and the com- 
panion's perturbing tidal gradient. The tidal param- 
eter /Z2 depends on the neutron star equation of state 
(EOS) and it is related to the relativistic generalization 
of the Newtonian, dimensionless, Love number [5j] k-i 
as fj,2 = 2/(3G)fc2i? 5 , where R is the star radius and 
G the Newton constant. The late-inspiral part of the 
GW signal, where tidal effects are stronger, can be used 
to measure the tidal Love number and thus to extract 
information about the nuclear EOS. A recent study @ 
of the measurability of G\xi, based on the tidal exten- 
sion of the effective-one-body (EOB) model [^.[8|-|Tl|. 
has shown that from a detection of GWs up to merger all 
normal matter content (npefi) EOS with maximum mass 
> 1.97 Mq can be distinguished at 95% confidence with 
signal-to-noise (SNR) ratio 16 and for any physical mass 
ratio El- 
Accurate theoretical modelling of GWs from BNS coa- 
lescence is a challenging task. High post-Newtonian (PN) 
tidal corrections and resummation techniques are needed 
to push the validity of the (semi) analytical models up 
to contact. Next-to-leading order (NLO, fractional 1PN 
accuracy) Q (then confirmed in [Hj]) and next-to-next- 
to- leading-order (NNLO, fractional 2PN accuracy) rela- 



On the contrary, if only the early inspiral waveform is considered, 
i.e. only GW frequencies < 450 Hz (for a BNS with individual 
masses 1.4 Mq) it is not possible to measure Gu2 with sufficient 
accuracy to discriminate among different EOS 0]. 



tivistic corrections to the tidal interaction energy have 
been computed recently using effective-field-theory tech- 
niques [l3| . Fractional 1PN tidal corrections to the wave- 
form were also obtained in [l4j]. The high-PN tidal 
corrections effectively amplify the magnitude of leading- 
order tidal effects, and are now incorporated in the tidal 
EOB model @, 0, H3 , which is currently the most sophis- 
ticated analytical tool available to model the dynamics 
and waveforms of neutron star (or even mixed) binaries 
up to contact. 

Numerical relativity (NR) simulations are the funda- 
mental tool to compute the dynamics and waveform of 
the last few orbits of a coalescing BNS system. NR data 
can be used to calibrate yet uncalculated higher-order 
tidal effects and to tests the reliability of the analytical 
models. To date, however, only few works have explored 
this important problem [16l4l8l | . 

A first comparison [l6| between waveforms from three- 
orbits NR simulations and the standard, point-mass, 
Taylor- T4 PN approximant, pointed out that the dephas- 
ing accumulated during the last orbits up to merger can 
be observed and used to constrain the EOS. As discussed 
there, a major limitation of that work was probably given 
by the length of the NR data available at the time. 

Long-term (nine and eleven orbits) BNS numerical 
simulations were presented in Refs. [l7], [Ha], and com- 
pared there with the prediction of the tidal EOB model. 
By performing a gauge-invariant and frequency-based 
analysis of the phasing, it was found that the tidal in- 
teraction predicted by the numerical simulation is im- 
portant even in the early part of the signal. To model it 
analytically, it was necessary to introduce effective frac- 
tional 2PN tidal corrections to the EOB model, yielding 
an amplification of the analytically predicted tidal ef- 
fects [7J, via the free parameter 0:2! the simulations con- 



2 



strained it0 in the range 40 < «2 < 100. That analysis 
also pointed out that a more detailed analysis ol finite- 
resolution uncertainties on long-inspiral BNS waveforms 
was needed to correctly estimate the magnitude of tidal 
effects. 

Recently, this important task was undertaken in [l8j . 
that presented the first comprehensive analysis of the 
uncertainties on the waveforms due to truncation and 
finite extraction error in a nine-orbit BNS simulation. 
In the same work, NR waveforms were compared to the 
tidal T4 app roximant [f| including also NLO tidal correc- 
tions (7l.ll2l.[l4l]. Significant effects on the phasing due to 
high-order tidal effects were observed in the "best data" 
(extrapolated from various resolutions), although a more 
conservative error estimate did not allow to distinguish 
higher-order tidal effects. 

Motivated by these works and by the last analytical 
results of we address in this paper the following 
question: Is the current analytical knowledge necessary 
and sufficient to reproduce NR data within their uncer- 
tainties? 

We present results about new NR long-term simula- 
tions and their comparison with the up-to-date EOB 
model u\. [Ta which includes the NNLO tidal correc- 
tions of [1J]. Two different sets of simulations of the 
same initial data are considered and some difficulties in 
obtaining NR data of sufficient accuracy are pointed out. 
We experimentally estimate the contact frequency of the 
binary, a fundamental information for a comparison with 
analytical models. NR data are then compared to the 
EOB model by carefully taking into account their uncer- 
tainties. We make use of the gauge-invariant relation 
between the (reduced) binding energy E and the (re- 
duced) angular momentum j of the system in order to 
analyse the dynamics of the binary [20|. The phasing of 
the waveforms is studied both in the time domain and 
by means of a gauge-invariant and frequency-based ap- 
proach which does not require to fix any relative (phase 
and time) alignment between the waveforms [l|| . 

The paper is organized as follows. In Sec. |TT]the tidal 
EOB model used here is reviewed. In Sec. Illll the numeri- 
cal simulations are presented. The EOB /NR comparison 
is discussed in Sec. lIVI Concluding remarks are in Sec.fVl 
We use units G = c = AI Q = 1, unless otherwise stated. 



II. TIDAL 2PN EOB MODEL 

The tidal extension of the EOB model of the binary 
dynamics has been defined in Ref. Q and then improved 
in Ref. [l3| and Q to include fractional 2PN corrections 
in the tidal part of the EOB potential A(r) . In particular, 
we address the reader to Appendix A of Ref. Q for a 



2 A similar conclusion, 012 ~ 40, was also reached in 0] using 
non-conformally-flat, NR stationary BNS sequencics. 



collection of ready-to-use formulas that define the EOB 
dynamics and waveforms including fractional 2PN tidal 
effects. Here we only summarize the main points. 
The EOB radial potential has the form 



A(u) = A u (u)+A tida \u), 



(1) 



where u = 1/r = GMj \c 2 r ab) is the Newtonian poten- 
tial, M = Ma + Mb the total mass, tab the relative 
separation, A°(u) denotes the point-mass potential and 
A tldal (u) is the supplementary tidal contribution of the 
form 
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The point-mass potential is defined using the usual Pade 
resummation of the 5PN Taylor expansion of the A func- 
tion with the 4PN and 5PN EOB parameters (05, a§), 
A°(u) = P^[l-2u + 2vu 3 + a i iyu 4: + a 5 vu 5 + a6iyu 6 }, where 
a 4 = 94/3-(41/32)tt 2 , v = M A M B /M 2 , and P™ denotes 
an (n,m) Pade approximant. Following the finding of 
Ref. Q (then substantially confirmed by Ref. [2ljj) we fix 
the free EOB parameters to the the values 05 = —6.37 
and 06 = +50 which lie within the extended region in 
the (05, ae) plane yielding a good fit of the binary black 
hole equal-mass simulations. 

In the tidal contribution, Eq. (|2|), the terms kJu 2£+2 
represent the leading-order (LO) tidal interaction, while 

the additional factor A[f^ al (u) takes into account the ef- 
fect of distance dependent, higher-order relativistic con- 
tributions to the dynamical tidal interactions: 1PN (first- 
order in u, or next-to- leading order, NLO), 2PN (of order 
u 2 , or next-to- next-to- leading order, NNLO), etc. The 
dimensionless EOB tidal parameter kJ is related to the 



tidal polarizability coefficients G[i 
star as 



A,B 



of each neutron 



(3) 



where 



(2£ - 1)!! 



M f 



M A (GM/c 2 ) 2 ^ 1 



(4) 



Here we take advantage of the new analytical results of 
Ref. [l3j and we use the newly computed expressions of 

^tidai( u ) for ^ = 2 and I = 3 at NNLO accuracy. Focus- 
ing on the most relevant equal- mass case (v = 1/4), the 
relativistic correction to the tidal potential reads 



A ( £l l (u) = l + a?>u + a 



(5) 



where the coefficients a^ 2 are i m this particular case, 
pure numbers. Specializing to the equal-mass case 
Eqs. (6.9)-(6.10) and Eqs. (6.21) and (6.22) of Ref. [3 
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we obtain 

i t ( Lw = i + !*+n u2 ' (6) 
A ( fLw = i+^+^ 2 : (7) 

43» - 1. (8) 

where we indicated explicitly the absence of (yet uncom- 
putcd) higher-order corrections to the t — 4 rclativistic 
contribution. The EOB waveform and radiation reaction 
is computed as in Refs. [l9| and takes explicitly into ac- 
count the 1PN tidal corrections of [l4[ (see Appendix A 
of [f| for the precise definition of the EOB waveform with 
tidal corrections). 

Equations <[6j)-((Sj) define the most advanced tidal EOB 
model based on analytical information only. In Ref. [l9| 
a slightly simplified representation of the functions A[ J j 
was used. Since at the time the NNLO calculation was 
not completed yet, and only a\ = 5/4 was known, one 
was using the following NNLO effective expression for the 
relativistic tidal corrections 

^Sal = 1 + "l« + «2« 2 , (9) 
(2) 

with a i = a\ = 5/4 fixed to be the same for I — 
2,3,4 and cti taken as a free effective parameter (for 
all I's) to be fitted for by comparison with NR simu- 
lations. Although in the following we shall mainly focus 
on the purely analytical 2PN tidal EOB model defined 
by Eqs. (EH)-©, we shall also contrast some predictions 
of the effective 2PN EOB model given by Eq. El with the 
numerical data. 

The last important concept we want to remember is the 
definition of contact between the two stars. This quan- 
tity is important because the analytical model ceases, in 
principle, to be valid after this moment. Such a formal 
contact moment was introduced in Eqs. (72) and (77) 
of 0] by the condition that the EOB radial separation 
tab becomes equal to the sum of the tidally deformed 
radii of the two stars, namely 

r co„tact = ^ + h^e A {r c lt^)) Ra + A^B, (10) 

where 6 a (j~ab ) = ^^BR A /( r AB^) 1S the dimensionlcss 
parameter controlling the LO strength of the tidal defor- 
mation of body A by its companion B, Ra is the star 
radius and hf is the shape Love number [l5|, [22[ • The 
dimensionless quantity is of order unity, but one is 
expecting it to be a function of the relative separation 
r, that increases as r decreases (this was found in a re- 
lated black- hole study ^3). In Ref. [l|| it was found 
that the effective value — 3 was necessary to allow 
this EOB-predicted contact to occurr always before the 
NR-defined merger. We shall briefly comment in Sec. IIVI 
about the magnitude of the amplification needed on hi 
so to reconcile the EOB contact defined by Eq. (fTTJ)) with 
a certain NR-defined contact. 



TABLE I. Summary of the grid configurations and of the 
runs. Columns: name of the configuration, maximum refine- 
ment level, minimum moving level, number of points per di- 
rection in the moving levels, resolution per direction in the 
level I — Z m ax, number of points per direction in the non- 
moving levels, resolution per direction in the level / = 0. 
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FIG. 1. (color online) Numerical quadrupolar gravitational 
waveform extracted at the outermost radius r D b s = 750 = 
247.55M for the CENO and WENO data (run H). Top: real 
part and amplitude (dashed lines). Bottom: frequency. The 
vertical lines mark the contact of the WENO data. 



III. NUMERICAL SIMULATIONS 

Target waveforms for the comparison with different an- 
alytical predictions are computed via NR simulations like 
those presented in great detail in [l8j]. While NR simu- 
lations of BNS have reached a certain maturity (see [23[ 
for the most recent review), intrinsic difficulties in the 
treatment of general rclativistic hydrodynamics (GRHD) 
make the numerical calculations of small effects in long 
simulations still very challenging. Note that the study 
of tidal effects in the late inspiral requires to resolve (at 
least) dephasing of < 0.5 rad over ten cycles. In par- 
ticular the numerical viscosity of high-resolution-shock- 
capturing scheme (HRSC) typically employed, is strongly 
dependent on the reconstruction scheme (cell-interfaces 
interpolation), and plays an important role in the accu- 
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FIG. 2. (color online) Numerical dynamics: E(j) curves 
for two series of simulation (CENO and WENO data), and 
for two point-mass analytical models (EOB resummed and 
"canonical" Taylor-expanded 3PN). The vertical dashed line 
marks the angular momentum at contact. 



racy of the simulations, see e.g. 0, Hf|. While partially 
under control in short (three orbits) runs by the use of 3rd 
order reconstructions, long term simulations such those 
presented here are challenging, also due to the computa- 
tional cost of extensive testing. 

In this work we consider long-term evolutions (~10 
orbits) of an equal-mass conformally-flat (CF) and ir- 
rotational initial configuration of Arnowitt-Deser-Misncr 
(ADM) mass Ma DM = 3.00506(2) and angular momen- 
tum Jadm = 9.716(1). The initial separation is d ~ 50 
associated with GW frequency of ~ 394 Hz. The fluid is 
described by a T-law EOS (r = 2) and isentropic evolu- 
tion were considered as in [18[ ■ The baryonic mass of each 
star is Mb = 1.62500(0), the gravitational mass of each 
star in isolation is M/2 = M A = M B = 1.51483(7), ra- 
dius and compactness are respectively R = 10.82065(0) 
and C = 0.14. The corresponding I = 2 dimensionlcss 
Love number is = 0.07890(1) and the £ = 2 shape Love 
number is h-2 = 0.8699. The initial configuration (2|| is 
computed with the LORENE library and publicly avail- 
able, and was already considered in Refs. |17l - fl9j . 

Evolutions were performed with the BAM code de- 
scribed in [25|, [13, ■ Here we mention the GRHD is 
handled with finitc-diffcrcncing HRSC based on primitive 
reconstruction, the Local-Lax-Friedrichs (LLF) central 
scheme for the numerical fluxes and Rungc-Kutta time 
integrators, see e.g. (2||. Cartesian grids and Berger- 
Oliger adaptive-mesh- refinement ("moving boxes" tech- 
nique) are used. The grid setup, resolutions, gauge pa- 
rameters, and finite differencing stencils for the metric 
sector are exactly the same as the convergent series dis- 
cussed in [l8|, they are listed for completeness in Table fl] 



We point out here that focusing on isentropic evo- 
lutions is justified by the following facts: (i) physi- 
cally, BNS evolutions are expected to be isentropic up 
to contact; (ii) any analytic model can not describe non- 
isentropic effects (e.g. shock heating) ; (iii) previous works 
demonstrated (l9j, |25j that considering non-isentropic 
evolutions actually leads to smaller tidal effects. 

Two series of simulations were performed: one is the 
convergent series presented in [lJI, where the HRSC 
employs the (formally) 3rd order convex-essentially-non- 
oscillatory (CEN03) for primitive reconstruction and a 
3rd order Rungc-Kutta scheme. The second series is 
computed with the same setup except the use of the 
(formally) 5th order weighted-essentially-non-oscillatory 
(WENOZ) method of [301 an d a Runge-Kutta scheme 
of 4th order. In [l8j it was presented a detailed anal- 
ysis of the uncertainties that affect the waveform due 
to truncation and finite-extraction errors; the new data 
computed for this work show analogous features. How- 
ever, we observed differences between the two data sets. 
For a given resolution, the merger in the CENO data 
occurs earlier than in the WENO data. The dominant 
multipole I = m = 2 of the metric waveforms, /122, ob- 
tained from the two different setup is displayed in Fig. [1] 
In abscissa we use the retarded time u = t — r*, where 
r* = r s + 2M ADM \og(r s / (2M ADM ) ~ 1), and r s is the 
Schwarzschild radius corresponding to the isotropic (co- 
ordinate) radius r. The waveforms are extracted at the 
outermost radius r = 750 = 247. 55M. The simulations 
compute waveforms from the Newman-Penrose scalar -04, 
that is then decomposed in spherical harmonics modes, 
ipl m . The metric multipoles hg m are calculated from the 
■0l m by integrating the relation ipf 11 — hi m . To do the in- 
tegration, we use the frequency-domain procedure of [3l| 
with a low-frequency cut off at ojq = 0.02/M. The signal 
is integrated from the very beginning of the simulation, in 
order to include also the initial burst of radiation related 
to the use of CF initial data. This radiation is often called 
(somehow improperly) "junk" radiation. Note that in the 
text for brevity we consider the metric waveform multi- 
plied by the extraction radius without explicitly chang- 
ing the notation, i.e. hi m = rhi m . As it is clear from the 
figure, at a given resolution, CENO and WENO wave- 
forms accumulate a significant relative dephasing towards 
merger; uncertainties due to the HRSC numerical viscos- 
ity become larger as the simulation time advances and 
eventually dominant over truncation (and finite extrac- 
tion) errors towards contact (Afo>22 ~ 0.07 [l8[, see be- 
low for an estimate of the GW frequency of the contact), 
where any convergent behavior is lost. For both data sets 
the higher the resolution, the later is the merger [l8|, [25| . 
In practical terms, Fig. [T] indicates that the differences in 
the HRSC effectively influence the magnitude of the tidal 
interaction between the two stars, from a larger value for 
the CENO data to a smaller one for the WENO data. 

This effect can be properly quantified by exploring the 
actual dynamics of the BNS system so to contrast it with 
the corresponding point-mass one. An intrinsic, gauge- 
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FIG. 3. (color online) Binary dynamics from WENO NR simulations. Contour plot of the rest-mass density on the equatorial 
plane. The snapshots indicate the contact happens around dynamical time t c ss 2382M. This dynamical time, corresponding 
to observer's retarded time u c = t — r° hB , locates the contact at GW frequency Mui22 ~ 0.078. Run H. 



invariant, way of doing so is by means of the relation be- 
tween the total energy £ and total an gula r momentum J 
of the system, £(J). Following Ref. |2(| that computed 
this quantity for binary black hole (BBH) systems, £ and 
J are obtained from the NR data as 



£ NH (u) = £0 



A cNR 
^rad 



(«), 



J NR (u) = \J -AJ™(u) 



NR; 



(11) 



(12) 



Here, £q = M^ DM /M and J = J^ DM /M 2 are the initial 
ADM mass and angular momentum expressed in units 
of the total gravitational mass of the stars in isolation; 
AfNR and AJ^ (expressed in the same units) are the 
radiated energy and angular momentum between the ini- 
tial (retarded) time u$ and u. They are computed from 
the multipole moments of the metric waveform hg m and 
of its time derivative hg m , as 



1 

16^ 



EE 

1=2 rn=0 Ju ° 

«„„ e 



du' 



he m (u') 



(13) 



A^ a R d («) = ^EE / du ' m ^ 



=2 m=l Ju ° 



hi m (u )h* em (u ) 



(14) 



In Eq. (|12[) we also included the x, and y component of 
the angular momentum loss, Aj7"N R and Aj7^ R , though, 
as expected, they are of order 10 -10 and negligible in 



practice. For convenience, we work with the binding en- 
ergy per reduced mass E = (£ — M)/fi and the dimen- 
sionless rescaled angular momentum j = J'/Mfx, where 
/j = M A M B /M. 

The main panel of Fig. [3] compares the NR relations 
E NR (j) computed for the two data series (CENO and 
WENO) with two analytical, point-mass, curves: the 
canonical PN expanded E(j) relation (see Eq. (5) of 01), 
dashed line, and the NR-tuned, EOB resummed one 
(dash-dotted line), that was found in [2(j to show an 
excellent agreement with corresponding BBH numerical 
curve (see Fig. 2 of 0, top panel). The E m (j) curves 
used here were obtained from waveforms taken at the 
the outermost extraction radius, r b s = 750 = 247. 55M, 
and, for simplicity, by including only the I = m = 2 mul- 
tipole in the calculation of Eqs. (fTTj) -(fT2" j) . To illustrate 
the influence of the uncertainties due to finite resolution, 
instead of displaying the NR data as a simple curve, we 
present them as the shaded band that is included between 
the medium (run M, bottom border) and high (run H, top 
border) resolutions. The diagram illustrates that, while 
the CENO data graze the EOB point-mass curve at the 
very beginning of the simulation (see inset) they visibly 
deviate from it after, indicating the presence of strong 



3 Although the canonical PN curve was very close to the CF 
initial state, it was found to progressively deviate from the nu- 
merical calculation, giving then an inaccurate representation of 
the point-mass dynamics. 
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tidal effects. On the contrary, the WENO data remain 
always very close to the point-mass EOB curve, so to be 
almost indistinguishable on the scale of the main plot. 

We experimentally conclude that the more dissipative 
numerical setup (CENO data) artificially amplifies tidal 
effects, leading to severe inaccuracies on the waveform 
phasing. The improvements obtained in the WENO data 
and the availability of convergence tests and error esti- 
mates [l8[ are crucial for the comparison with the ana- 
lytical information presented below. 

Let us finally briefly comment the inset of Fig. [2j which 
focuses only on the initial part of the i? NR (j) curve. Sim- 
ilarly to the black- hole case (see inset of Fig. 1 of 20]), 
the initial state of the system is very close to the point- 
mass, 3PN canonical E(j) curve. Then the effect of the 
losses due to the junk radiation moves the initial state 
down, close to the EOB curve. Note in addition that 
this early-inspiral part the NR curve is above the point- 
mass EOB curve because only the I = m = 2 mode has 
been included. We shall explore the effect of the other 
multipolcs in the next Section, in the context of the de- 
tailed comparison with the tidal EOB model. 

For the comparison with analytical predictions, it is 
important to have an estimate of some contact frequency 
extracted from the NR data, since the tidal EOB model 
ceases to be valid after. However, connecting the (local, 
strong-field) dynamics, parametrized by the dynamical 
time t, of the two objects with the radiation observed in 
the wave-zone unambigously, is a non-trivial task. In first 
approximation, any phenomenon occurring in the strong 
field region at dynamical time t reaches the observer r b 8 
at retarded time u = t — r* (r bs)- For simplicity, in the 
following we shall assume such relation to connect the 
two events, though, in doing so, we are neglecting an ad- 
ditional time-delay due to the propagation of the signal 
in the strong-field region. In Fig. [3] we show snapshots of 
the rest-mass density in the orbital plane at few dynam- 
ical times around the contact of the two stars. The pic- 
tures indicates clearly that at t c « 238271/ the well-known 
shearing contact is taking place. The corresponding GW 
frequency is Mu 2 2{u c = 2382M) « 0.078, and the cor- 
responding value of the angular momentum is j c ~ 3.63. 
The latter is shown in the dashed vertical line in Fig. m 
Note that the final merger (e.g, formally defined by the 
first peak of I/122I) occurs later, at Muj 22 ~ 0.13. 



IV. EOB/NR COMPARISON 

In this Section we consider only the WENO data and 
compare them with the tidal EOB model. We present 
two types of comparisons, one for the dynamics, through 
the E(j) curve, and one for the phasing. We shall take 
as "best" (multipolar) waveform the one computed with 
the highest resolution available and extracted at the out- 
ermost radius r D b s = 247. 55M. 



A. Dynamics 

Let us first compare the NR E NR (j) relation to the 
corresponding analytical prediction. This is done in the 
left panel of Fig. [4] Together with the point-mass ana- 
lytical curves (EOB and canonical 3PN) already shown 
in Fig. |31 we also display the analytical E(j) EOB curve 
that includes tidal effects (solid, red online) at NNLO 
order. The numerical curve (solid, black online) is dif- 
ferent from the corresponding band of Fig. [3] in that 
and A£ r ^ R include all multipolar contributions 
up to f max — 4. Note that, although the symmetry of 
the system implies that all multipolcs with odd values of 
m have to be zero, they are actually nonzero, and are 
included in the computation of the losses to which they 
contribute as very small amplitude, structureless, noise. 
In these simulations, differently from (l9| . we do not im- 
pose rotational ('V) symmetry on the orbital plane, but 
evolve instead the equations in the whole plane z > 0. 
To give numbers in the t = 2 case, during the inspiral it 
is, at maximum, | /z-22 1 ~ 5 x 10~ 3 , \h 2 o\ ~ 2 x 10 -5 and 
\h 21 \ - 6 x 10~ 6 . 

Though the m = modes are practically negligible 
during the inspiral, they actually contribute to the ini- 
tial burst of radiation and thus they must be included 
in the computation of E NR (j). For example, during this 
epoch, that lasts for <~ 25M, | /120 1 ~ 0.5 x | ] and 
I A40 1 ~0.5x|/i44|. By contrasting the insets of Fig. [3] and 
of the left-panel of Fig. |H we observe that the subdom- 
inat multipoles drive the NR curve not only below the 
point-mass curve, but even below the tidal EOB curve El- 
Like in Fig. |31 the dashed vertical line marks the location 
of the NR-contact. We note in passing that to make the 
EOB-contact consistent with the NR-defincd contact we 
should replace in Eq. ([TUf the value h 2 = 0.8699 (which 
would give a EOB contact frequency slightly smaller than 
the merger 

frequency Mu 22 ~ 0.094) with hf ~ 3.4, 
which yields an EOB contact frequency ~ 0.078. This es- 
timate of the deformation of each star coming from actual 
data is consistent with the guess of Ref. [19[ {h 2 S = 3), 
but further work will be needed to understand, numeri- 
cally, the actual amplification experienced by h 2 . 

By close inspection of the plot we observe that the NR 
curve actually change its slope: around j ~ 3.8 it lies 
between the tidal and point-mass curves, while during 
contact the curve bends again below the tidal EOB one. 
We argue this behavior is caused by numerical inaccura- 
cies, for example small violations of mass conservation, 
to which the computation of £' NR (j) is extremely sensi- 
tive. It was not possible to identify the precise cause, but 



4 The tidal EOB curve is below the corresponding point-mass one 
indicating that, due to the attractive nature of tidal interaction, 
the system is gravitationally more bound. 

5 Similarly, the "bare" contact, with /12 = 0, is even closer to the 
merger, giving M1U22 ~ 0.106. 
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FIG. 4. Comparison between EOB and NR dynamics. Left panel: reduced binding energy (E) versus reduced angular 
momentum (j) curves. Right panel: differences with the tidal EOB model. The shaded region represents the estimated 
uncertainties on the "best" numerical curve. Numerical data are consistent with the tidal EOB model. 



wc mention that it corresponds to a numerical oscillation 
at the fourth digit of £ NR (j). 

The differences between numerical data and the an- 
alytical predictions are made precise in the right panel 
of Fig. |U The four curves represent the four differences 
E x - E EOBtidM where the label X indicates in turn: NR 
(thick, red online) with £ max = 4; NR (thin, black on- 
line) with £ max = 2; NR (thick, dashed, black online) 
with £ max = 4 extrapolated at infinite extraction radius; 
EOB point- mass (dash-dotted, black online). The shaded 
region represents the estimated uncertainty affecting the 
best waveform. It has been obtained by taking into ac- 
count the following three effects: (i) resolution; (ii) finite 
extraction radius; (iii) contribution of higher multipolcs. 
To compute it, we first took the following differences be- 
tween E NR (j) curves: (i) between H and M data, AE HM ; 
(ii) between H data extrapolated at infinite extraction ra- 
dius and at r b s = 247.55M, AE^. The extrapolation 
was done from data at r ohs = {400, 500, 600, 700, 750}, 
aligned using retarded time, taking a quadratic fit; (iii) 
between H data with £ max = 4 and with £ max = 2, AEf. 
Then the final (conservative) error-bar is obtained as 

A NR E(j) = ±^{AEgf + (A£ HM ) 2 + {AEff. Note 
that at large values of j (beginning of the simulation) the 
most relevant uncertainty is the one due to the finite ra- 
dius and to the choice of £ max ; after that, finite-resolution 
effects become dominant. 

The analysis of the dynamics shows that the NR data 
are consistent with the state-of-the-art tidal EOB model, 
and distinguishable from the point-mass EOB model up 
to contact (j c = 3.62). Remarkably the NR curves stay 
very close to the tidal EOB up to j ~ 3.5, point at which 
the differences with the point-mass curves are at least 
two-sigma beyond the uncertainties. 



TABLE II. Phase differences A0 = EOB - NR , with its 
uncertainty cta<p, accumulated between frequencies [wijWg] = 
[0.041, 0.062] /M obtained by integrating the difference be- 
tween the NR Qui and some of the EOB curves of Fig. [5] 
the point-mass, the tidal LO and the 2PN (NNLO) analyti- 
cal one with c4 2) = 85/14. The error-bar ±<ta0 is obtained 
by integrating the shaded region around the NR Q u curve in 
Fig. © 



EOB model 


point-mass 


LO tidal 


NNLO tidal 


A(/> [rad] 


3.92 


1.49 


1.06 


cta<p [rad] 


0.61 


0.61 


0.61 



B. Waveforms and phasing 

Let us now compare the waveforms and quantify the 
phase difference accumulated between NR and EOB data 
in the last orbits of the inspiral. Focusing on the £ = m = 
2 mode, we first present a time domain comparison of 
the amplitude and frequency, and then switch to a more 
quantitative analysis of the phasing. Higher multipoles 
are discussed at the end of the section. 

In Fig. [5] we contrast the NR waveform (thin solid line, 
black online) with 4 different analytical models: the EOB 
point-mass (no tides, dash-dotted line), the tidal EOB 
with only LO tidal effects (i.e., we set A" dal = 1 for all t s 
in Eq. (|2J|), dashed, thick curves, blue online), the 2PN- 
accurate (NNLO) tidal EOB model (solid thick curves, 
red online) and the PN-cxpanded Taylor T4 model with 
leading order tidal effects (see Sec. IIIC of [19| ). We show 
together the gravitational frequency (top left) and wave- 
form modulus (top right) while the corresponding dif- 
ferences (EOB-NR) are exhibited in the bottom panels. 
The vertical dashed line on each panel indicates the NR 
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u/M u/M 





u/M 



u/M 



FIG. 5. Comparing numerical and analytical I = m = 2 waveform: frequency (left panels) and modulus (right panels). The 
vertical line locates the NR contact. The shaded regions in the bottom panels are error estimates on the NR frequency and 
modulus. 



contact frequency Mu^ ~ 0.075, after which we don't 
expect any analytical model to be accurate. The relative 
time, t and phase, a, shifts that are necessary to align 
the analytical to the numerical waveform are determined 
using the procedure described in Sec. VA of Ref. filil ]. 
that relies on the minimization of the \ 2 °f the phase 
difference over a certain frequency interval. Here we use 
the frequency interval (cji,lj 2 ) = (0.038, 0.049)/Af, (cor- 
responding to (u 1 ,U2)/M = (165.8,1610.9), that begins 
after the initial burst of radiation, so as to remove also 
possible inaccuracies due to the integration procedure 
needed to get /122 from ^| 2 . The shaded regions in the 
bottom panels indicate the NR uncertainty, that takes 
into account finite-extraction radius and finite resolution 
effects. It is obtained by: (i) Richardson extrapolating 
the two highest resolutions assuming second order con- 
vergence and taking the difference with run H data; (ii) 
similarly, taking the difference between the NR wave- 



form extrapolated in extraction radius H and the one at 
^obs = 750. In practice, for the first half of the simula- 
tion (u ~ 1000M) the uncertainty is dominated by finite- 
extraction-radius effects, while later it is the resolution to 
play the most important role. The two contributions are 
then summed in quadrature and one takes the (±1/2) of 
the square root to obtain a two-sided error bar. Note that 
the finite-resolution uncertainties we quote are consistent 
with Table II of [l8| and provide an average between 
optimistic and conservative estimate of the errors (ob- 
tained respectively by resolution-extrapolated data from 
five and three simulations with different resolutions). 

From the top panel of Fig. [5] one sees visually how 
the tidal models clearly yield a better agreement with 



6 In Ref. Il8l it was shown that in this case phase difference due 
to finite-radius effects varies between ~ 0.2 rad during the early 
inspiral to ~ 0.1 rad at merger, while the fractional difference in 
amplitude varies between 1% and 0.5%. 
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-o-WENO (H) 




WENO uncertainty 




EOB point-mass 




EOB tidal (LO) 




EOB tidal 2PN af ] = 85/14- 




EOB tidal 2PN eir a 2 = 20 




EOB tidal 2PN L,fr a 2 = 40 




EOB tidal 2PN'' ,r a 2 = 100 







FIG. 6. (color online). Comparison between various EOB 
Qu; curves and the NR one. The good visual agreement be- 
tween the LO and 2PN tidal EOB models and the NR curve 
highligths the presence of tidal interaction during the inspi- 
ral. Note that the effective value «2 = 100 used in Ref. [lQI 
is incompatible with the NR curve. 



0.(14 
0.035 
0.03 





0.04 
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0.03 
^ 0.025 
3 0.02 
0.015 
0.01 
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WENO u 
-•EOB (point-mass) 
-EOB ridid LO 
-EOB tidal 2PN a%> = 85/14 
-EOB tidal 2PN* B a a = 20 
■ EOB tidal 2PV fl ai = 40 
- T4 tidal LO 



u/M 



FIG. 7. (color online) Phase difference between various an- 
alytical phase and the NR phase. The vertical line indicates 
the NR contact. The phase difference of ~ — 0.2rad between 
the 2PN tidal EOB model and the NR waveform is compatible 
with the error estimates on the latter. 



the numerical data than the simple point-mass EOB 
model Q- This information is made more quantitative in 
the the bottom panels, where the point-mass analytical 



Note that up to GW frequency M0J22 ~ 0.1 the Taylor T4 point- 
mass phasing agrees very well with the NR one and thus with 
the NR-tuned EOB point-mass. 



FIG. 8. (color online). Comparing EOB and NR waveform 
moduli for the most relevant subdominant multipoles /123 and 
/144. The vertical line marks the NR-defined contact time. 



prediction is seen to deviate away starting from 2000 M, 
while the other differences remaining much flatter, and 
marginally close to the error bar, up to contact. The 
comparison of Fig. [5] allows us to deduce the presence 
of tidal effects in the very late part of the inspiral, just 
before contact. However, it also indicates that this com- 
parison is sensible at most to leading-order tidal effects, 
since, given the uncertainties on the NR data, it is not 
possible to meaningfully disentangle 1PN and 2PN tidal 
corrections from the leading-order ones. On these plots, 
the Taylor T4 PN model with LO tidal corrections looks 
consistent with the EOB predictions, yielding similar dif- 
ferences with the NR waveform. Note that at the very 
contact position the frequency and amplitude differences 
show a clear increasing trend. This might either be due to 
the "blurred" nature of contact in the NR data, or to the 
lack of suitably determined next-to-quasi-circular correc- 
tions. For simplicity, we will not overtune here more our 
analytical model and postpone to a future investigation 
the detailed analysis of these additional effects. 

We discuss now the phasing by means of a gauge- 
invariant and frequency-based analysis employing the 
Q u = u! 2 /lj function (where we put lu = Mu>22 for 
simplicity) introduced and used extensively in Ref. [l9j . 
We recall that the meaning of this function is that the 
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time-domain GW phase 0( Wl . W2 ) accumulated between 
frequencies (wi,^) is given by the integral 4>( u , 1 ^ 2 ) = 
Qujdlnu). Consequently, a change of Q u (uj) of the or- 
der ±1 during a frequency octave ln^/wi) corresponds 
to a local dephasing (around ui) of 8(f) ~ ±1 rad. The 
main advantage of the diagnostics is that it is in- 
dependent of the arbitrary time and phase shifts (r, a) 
necessary to compare the waveforms in the time domain. 

As in Ref . [l9[ we cannot compute from the raw NR 
data, but we have fitted the phase of /122 with a suitable 
PN-based expression (see Eqs. (27)-(28) of pj). Here 
the best fit is given by using a sixth-order polynomial (in 
contrast to the fourth-order polynomial employed in (l9l |) 

in the variable x = [v{t m — t)/5] (where t m is a fit- 
ting parameter formally representing the merger time), 
and the time interval [t^, tji]/M = [965, 2400], which cor- 
responds to frequencies [ujl,wr] = [0.042, 0.063]/A/. The 
curve, so obtained, is represented by the thick-solid 
line with circles in Fig. [51 The shaded region around 
the curve indicates the uncertainty on the curve as given 
by i ctqnr. This numerical uncertainty Oqnr is 

estimated by putting together the effect of truncation 
error, of finite extraction radius, and of the fit. To do 
so, we first calculated other three curves,: one from 
the waveform of the M run, Q^ Rm ; another from the 
waveform of the H run, but extrapolated at infinite ex- 
traction radius, Q^ Rh ; a third one doing the fit of the 
H data with the fourth-order polynomial instead of the 
sixth-order one, Q w RH,n=4 . We then computed the differ- 
enecs SQu'°° = QT H - Qu R ", 6Q™ = Q™» - Q* Rm 
and <5Q2'™= 4 = Q£ Rh - Q^ Rh '" =4 and summed them 
in quadrature, so to estimate the error-bar ctqnr = 

±1/2^ (6Q™ 1 ) 2 + (6QK°°f + ((5gS'" =4 ) 2 that we rep- 

resented in Fig. [5] This conservative error estimate is of 
order unity, as it varies between (—3, +2) at io = 0.0415, 
and between ( — 1, +1) at ui = 0.0603. 

Together with the numerical curve we also exhibit in 
the picture other five analytical EOB curves: the point- 
mass EOB (dash-dotted, black online), the LO tidal EOB 
(thick-dashed, blue online), the analytical, 2PN tidal 

EOB with a 2 ' =85/14 (thick-solid, red online) and the 
effective 2PN tidal EOB with the effective values a.2 = 40 
and OC2 = 100 used in Ref. [l9| (lowermost dotted lines, 
black and magenta online, from top to bottom). The 
figure highlights clearly the dependence of NR on tidal 
interaction. One concludes that: (i) the NR curve is al- 
ways very close to the 2PN tidal EOB one; (ii) it is very 
well distinguishable from the point-mass prediction; (iii) 
the error-bar on the NR curve is too large to appreci- 
ate the differences between the LO and 2PN tidal EOB 
models; (iv) the effective 2PN tidal EOB model used in 
Ref. [l9| with a.2 = 100 significantly overestimates the 
magnitude of tidal interactions in the NR data; (v) the 
effective 2PN tidal EOB model with a.2 = 40 gives a good 
average of the numerical points. 



The most important information suggested by Fig.|5]is 
that the EOB tidal model constructed using only analyt- 
ically computed tidal information is by itself consistent 
with the NR simulation, without the need of tuning any 
additional tidal EOB flexibility parameter yielding an ef- 
fective amplification of the tidal interaction as the stars 
get closer and closer. We can not exclude that such an 
amplification exists 0, but one will need much higher ac- 
curacy in the late inspiral phase to identify actual phys- 
ical effects. In this respect, although the value a.2 = 40 
fits well the NR data, it does not indicate definitely an 
amplification of tidal effects, as truncation errors are still 
dominant in this frequency range and the value of a.2 is 
very sensitive to small changes on that are barely vis- 
ible on the plot. For instance, note how the EOB is 
easily matching the upper bound of the error-bar by tak- 
ing a.2 = 20. We expect that the use of higher- resolutions 
and/or more accurate numerical treatments of the hydro- 
dynamics will further move up the NR Q u curve, so to 
favour smaller values of a.2 than larger ones. 

The (visually) small differences between the Q^'s in 
Fig. [6] actually correspond to relevant dephasings, of or- 
der 1 rad or more. This information, relative to the fre- 
quency interval [uji,uj2]M = [0.041,0.062], is quantified 
in Table HH The 2PN (NNLO) tidal EOB model accu- 
mulates a dephasing of 1.06 rad, the LO EOB model 
1.49 rad, while the point-mass EOB 3.92 rad. The un- 
certainties on these numbers are of the order £ta</> = 
0.61 rad, and obtained by integrating the shaded region 
in Fig. El 

Due to the fitting procedure involved in the computa- 
tion of the curves [l9j . it is important to verify the 
phasing with another diagnostic. Hence, we present also 
an analysis based on waveform alignment as customary in 
the literature. The time evolution of the phase difference 
A(t> 22 (t) = (j) x (t) - NR (i) (where the label X can be ei- 
ther EOB or T4) is shown in Fig. [7] It is computed from 
the time-and-phase alignement waveforms, as in Fig. [5] 
The shaded region represents the uncertainty on the NR 
phase. The qualitative information given by this plot 
confirms the analysis of the phasing given by the Q^: 
tidal effects are clearly visible before contact and the cur- 
rent analytical knowledge is sufficient to match the NR 
phasing up to contact (dashed vertical line in the figure) . 
It is not possible, however, to distinguish in the NR data 
the effect of higher-order tidal effects from LO ones. In- 
terestingly, on this plot the T4 tidal LO model performs 
marginally worse than the corresponding EOB model, 
~ —0.3 rad at contact and notably out of the estimated 
NR uncertainty. Note that the phase difference varies in 



Especially on the basis of the analytical considerations of 
Ref. [iii suggesting that might be replaced by an effec- 

tive distance dependent coefficient a^{u) = off /(I — ^lr^Oj 
tlr denoting the EOB effective light-ring location. Note how- 
ever that the corrisponding Q u curve would be indistingushable 
from the nonrcsummcd one on the plot. 
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the range ±0.2 rad from the beginning of the simulation 
(after the initial burst) up to contact. This value is con- 
sistent with lower value A0 — cta</> ~ 0.5 obtained from 
Table [II] For completeness, in Fig. [7J we also added the 
two phase differences with the effective 2PN tidal EOB 
model, with a.2 = 20 and a.2 = 40. The are both well 
within the error bar, although, consistently with the Q u 
analysis, the value a.2 = 20 actually yields a smaller de- 
phasing at contact. In conclusion, putting together the 
information of Figs. [1)][7J state-of-the art numerical simu- 
lations allow us to conclude that, if any actual amplifica- 
tion of tidal effects exists, it yields, conservatively, a < 40 
(as a conservative estimate), or, more likely, a>2 < 20. 

We conclude this Section with a few comments about 
the accuracy of the higher-order multipoles, that are ac- 
tually included in the computation of E NR (j). Figure [5] 
compares the modulus of the most relevant subdominant 
numerical multipoles £ — 3, m — 2 and £ = m = 4 with 
the corresponding EOB waveforms. To our knowledge, 
this comparison has never been shown before. The vi- 
sual agreement is rather good, with both analytical mul- 
tipoles averaging the corresponding numerical ones prac- 
tically up to contact, as it was the case for the £ = m = 2 
case. From the picture one also sees the large initial burst 
of junk radiation that must be included in the accurate 
computation of the E Nn (j) relation. 



V. CONCLUSIONS 

In this paper we have presented a comparison between 
dynamics and waveform from BNS coalescence computed 
from long-term (~ ten orbits) NR simulations and the 
tidal EOB model including all the known tidal PN cor- 
rections [l3j |. 

New numerical simulations have been presented which 
improve quantitatively previous results [18| . A set of sim- 
ulations which em ploy the same initial data, grid setup, 
and resolutions of [18j . but adopt an higher-order recon- 
struction method in the HRSC scheme and time integra- 
tor, has shown that tidal effects can be overestimated 
by numerical inaccuracies. While the data show conver- 
gent behavior before the contact, around this point and 
later the uncertainties related to numerical viscosity do 
not seem completely under control, and actually become 
dominant over truncation errors. 

In order to compare NR data with the analytical EOB 
model, we estimated the GW frequency of the contact 
as, Mu c 22 ~ 0.078. 

The dynamics of the system has been investigated by 
means of the E(j) relation between the reduced binding 
energy E and the reduced angular momentum j, com- 
puted here for the first time for BNS simulations and 
presented in Fig. |U The tidal EOB model is consistent 
with the NR data up to contact (j = 3.63) and even later, 
up to j ~ 3.5. 

The effects of tidal interactions are clearly visible in 
the NR/EOB waveforms. The comparison of amplitude 



and frequency (Fig. [5]) indicates that tidal effect become 
significant in the very late part of the inspiral, just be- 
fore contact. Given the uncertainties on the NR data, it 
is not possible to meaningfully disentangle 1PN (NLO) 
and 2PN (NNLO) tidal corrections from the LO ones. 
The T4 tidal model with LO tidal corrections is slightly 
worse (-0.3 rad) than the EOB model when getting close 
to contact. The phasing was studied by means of both a 
gauge- invariant, frequency-based analysis employing the 
diagnostic |19j and a standard time and phase-shift 
alignment procedure. The results of the former method 
are collected in Fig. [5] and Table HTl the ones of the latter 
in Fig. [7] The diagnostic is more affected by the noise 
of the data which results in somehow larger uncertainties; 
the time and phase alignment suffers of ambiguities in the 
choice of the interval and potentially underestimates ac- 
tual differences accumulated up to the alignment interval. 
In summary, while they give slightly different numbers, 
the picture emerging from the two phasing analyses is 
consistent: tidal effects are clearly visible during the late 
inspiral up to contact and the current analytical knowl- 
edge is sufficient to match the NR phasing. It is not 
possible, however, to distinguish in the NR data the ef- 
fect of higher-order tidal effects from LO ones. We ob- 
serve that after contact and up to merger (i.e., for about 
one further orbit), the 2PN EOB model performs better 
than any other analytical tidal model; note however that 
the extension of any analytical model beyond contact has 
only an effective meaning. 

In conclusion, we have shown that the current analytic 
knowledge incorporated into the EOB model is sufficient 
to reproduce within the uncertainties the numerical data 
up to contact. No calibration of any tidal effective-one- 
body free parameter is required, beside those already fit- 
ted to binary black holes data. While the 2PN (NNLO) 
model minimizes the differences with the NR data, it is 
not possible to significantly distinguish it from the 1PN 
(LO) model. Obviously, we cannot exclude the presence 
of a further amplification of tidal interaction as the star 
gets close (as suggested by Ref . [l3| ) , but the present NR 
data indicate that this effect, if present, is smaller than 
what believed in the past [l!| (i.e. a 2 = 100) and it is not 
possible to estimate it precisely. A conservative analysis 
points points to a.2 < 40, though we think that a more 
likely estimate (at one-sigma level) is ct2 < 20. Note 
in addition that for higher, more realistic compactnesses 
(say C = 0.16 — 0.18) tidal effects are even smaller, thus 
potentially more difficult to extract from the numerical 
data. Similar considerations also hold for the use of real- 
istic EOS, which present their own numerical challanges 
to be used in NR simulations. 

As a consequence, the 2PN-accurate tidal EOB 
model [HI used in this work@ should be considered in the 
future as the most reliable choice to produce exact/target 



Here we focused on the equal-mass case only, but we expect 
this EOB model to be accurate in the range of plausible mass 
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data for the development of templates for data-analysis 
purposes @. 

This work also pointed out the importance of exten- 
sive numerical tests to assess the uncertainties of the nu- 
merical data, and the potential need of new numerical 
strategies to perform accurate simulations. Considering 
that the simulations presented here are the longest and 
employ the among the highest resolutions to date, error 
assessment and convergence tests appear absolutely nec- 
essary in future studies of this kind. Because the use of 
significantly higher resolutions (e.g. ~ 400 3 points cov- 
ering each star) and extensive tests seem to be compu- 
tationally unfeasible, the development of alternative and 
more accurate numerical methods seems unvoidable to 



further improve and confirm our results. 
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ratios for a BNS system, from Ma/M b = 1 to Ma/M b = 0.7, 
corresponding to u £ [0.2422, 0.25]. The reason for this being the 
fact that the EOB model used here was found consistent with 



(now relatively old) BBH data up to 2:1 mass ratio (y = 2/9 = 
0.2222). Further improvements of the current EOB point-mass 
model are currently in progress |32| , consistently with the finding 
of Rcf. [H. 
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